## REPLICATION CODE FOR: Particularised Preferences for Civilian Protection
# Sophia Hatz and Lisa Hultman
# Code by Sophia Hatz

# Replicates Figures 1--6 in the manuscript

# Set up ------------------------------------------------------------------
set.seed(1237658)
rm(list = ls())


source("PP_Replication_Packages_and_Functions.R")

load("PP.RData")
load("PP_Labels.RData")


# global theme options for plots
theme_set(theme_bw() +
theme(
plot.title = element_text(size=10,face = "bold"),
axis.title.x = element_text(size=10,face="bold"),
axis.text.x = element_text(size=8,face="bold"),
axis.title.y = element_text(size=10),
axis.ticks = element_blank(),
legend.position = "bottom", 
legend.title = element_text(size=10,face = "bold"),
legend.background = element_rect(colour = "black"),
strip.background = element_blank(),
strip.text = element_text(face="bold", size=10) ))




# Figure 1: Support for civilian protection across six proposals ----------

outvars6 <- c("Supp_Aid","Supp_Obs","Supp_Pol","Supp_Mil","Supp_Ref","Supp_Rep")
outvars6.lab <- c("Humanitarian aid", "UN observers","UN police","UN military troops","Accept refugees","Repatriation")

PP$num <- 1:nrow(PP)
df <- reshape2::melt(PP[,c("num",outvars6)],
                     id.vars = "num",
                     measure.vars = outvars6,
                     variable.name = "variable") 
dfwc <- summarySEwithin(df, measurevar="value", withinvars="variable",  idvar="num", na.rm=TRUE, conf.interval=.95) 
dfwc <- dfwc %>% mutate(variable = factor(variable,labels=paste(outvars6.lab,"\nn= ",dfwc$N)))
df <- df %>% mutate(variable = factor(variable,labels=paste(outvars6.lab,"\nn= ",dfwc$N)))

png(filename = "Figure_1.png")

ggplot(data= df, aes(y= value, x= variable)) + 
  geom_point(data= df, aes(y=value),
             position = position_jitter(width = 0.2, height = 0.1),
             alpha = 0.01) +
  geom_errorbar(data= dfwc, aes(ymin = (value - ci), ymax = (value + ci)), width = 0.05) + 
  geom_point(data= dfwc, aes(y= value, x= variable), shape=21, size=1, fill="white") +
  xlab("") + ylab("Level of support") +
  geom_text(data=dfwc,aes(x= variable,y=value,label=paste("mean=",round(value,2))),size=3,nudge_y=-0.5) 

dev.off()



# Figure 2: Effect of identity factors on support for civilian protection --------

## define outcome variables used in hypothesis tests: five main outcomes plus index
(outvars <- outvars7[c(1:5,7)])
(outvars.lab <- outvars7.lab[c(1:5,7)])

## run models for each outcome variable and each identity factor to get p-values to adjust 
plist1 <- list()
plist2 <- list()
plist3 <- list()

# identity factor 1
p.values <- NA
estimates <- NA
for (i in 1:length(outvars)){
df <-
    PP[,c(indvars[1],outvars[i])] %>%
    gather(dv, value, outvars[i]) %>%
    mutate(IDV=factor(get(indvars[1]),labels=c("OW","Swedish citizens")))
  
df <- broom::tidy(estimatr::lm_robust(update(covadj_model, get(outvars[i]) ~ get(indvars[1]) + .), data = PP)) %>% filter(term=="get(indvars[1])") %>% mutate(IDV = indvars.lab[1], DV = outvars.lab[i]) 
  estimates[i] <- df$estimate # store estimates 
  p.values[i] <- df$p.value # store p-values
  
plist1[[i]] <- ggplot(df, aes(x = IDV, y = estimate)) + facet_wrap( ~ DV, scales = "fixed", shrink = TRUE) + geom_point(size = 2) +
    scale_y_continuous(limits=c(-0.28,0.25)) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
    geom_hline(yintercept = 0, colour = gray(0.5), lty = 2) + ylab("") + xlab("") + theme(axis.text.x = element_blank())
}

# add adjusted p-values
p.adj <- p.adjust(p.values[1:5], method = "holm") # N = total number of tests (DVs) for hypothesis 1
p.adj <- c(p.adj,p.values[6]) # use unadjusted p-value for index

for (i in 1:length(outvars)){
  df <- data.frame(p=p.adj[i],estimate=estimates[i])
  plist1[[i]] <- plist1[[i]] +
    geom_text(data=df,aes(x=1,y=estimate,label=paste("p=",round(p,2))),col="black",size=3,nudge_y=0,nudge_x = 0.3) +
    geom_text(data=df,aes(x=1,y=estimate,label=paste("Diff=",round(estimate,2))),col="black",size=3,nudge_y=0.05,nudge_x = 0.3)
}

# identity factor 2
p.values <- NA
estimates <- NA
for (i in 1:length(outvars)){
  df <-
    PP[,c(indvars[2],outvars[i])] %>%
    gather(dv, value, outvars[i]) %>%
    mutate(IDV=factor(get(indvars[2]),labels=c("OW","Women and children")))
  
  df <- broom::tidy(estimatr::lm_robust(update(covadj_model, get(outvars[i]) ~ get(indvars[2]) + .), data = PP)) %>% filter(term=="get(indvars[2])") %>% mutate(IDV = indvars.lab[2], DV = outvars.lab[i]) 
  estimates[i] <- df$estimate # store estimates 
  p.values[i] <- df$p.value # store p-values
  
  plist2[[i]] <- ggplot(df, aes(x = IDV, y = estimate)) + facet_wrap( ~ DV, scales = "fixed", shrink = TRUE) + geom_point(size = 2) +
    scale_y_continuous(limits=c(-0.28,0.25)) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
    geom_hline(yintercept = 0, colour = gray(0.5), lty = 2) + ylab("") + xlab("") + theme(axis.text.x = element_blank())
}

# add adjusted p-values
p.adj <- p.adjust(p.values[1:5], method = "holm") # N = total number of tests (DVs) for hypothesis 1
p.adj <- c(p.adj,p.values[6]) # use unadjusted p-value for index

for (i in 1:length(outvars)){
  df <- data.frame(p=p.adj[i],estimate=estimates[i])
  plist2[[i]] <- plist2[[i]] +
    geom_text(data=df,aes(x=1,y=estimate,label=paste("p=",round(p,2))),col="black",size=3,nudge_y=0,nudge_x = 0.3) +
    geom_text(data=df,aes(x=1,y=estimate,label=paste("Diff=",round(estimate,2))),col="black",size=3,nudge_y=0.05,nudge_x = 0.3)
}

# idv 3
p.values <- NA
estimates <- NA
for (i in 1:length(outvars)){
  df <-
    PP[,c(indvars[2],outvars[i])] %>%
    gather(dv, value, outvars[i]) %>%
    mutate(IDV=factor(get(indvars[2]),labels=c("OW","IS-affiliation")))
  
  df <- broom::tidy(estimatr::lm_robust(update(covadj_model, get(outvars[i]) ~ get(indvars[3]) + .), data = PP)) %>% filter(term=="get(indvars[3])") %>% mutate(IDV = indvars.lab[3], DV = outvars.lab[i]) 
  estimates[i] <- df$estimate # store estimates 
  p.values[i] <- df$p.value # store p-values
  
  plist3[[i]] <- ggplot(df, aes(x = IDV, y = estimate)) + facet_wrap( ~ DV, scales = "fixed", shrink = TRUE) + geom_point(size = 2) +
    scale_y_continuous(limits=c(-0.28,0.25)) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
    geom_hline(yintercept = 0, colour = gray(0.5), lty = 2) + ylab("") + xlab("") + theme(axis.text.x = element_blank())
}

# add adjusted p-values
p.adj <- p.adjust(p.values[1:5], method = "holm") # N = total number of tests (DVs) for hypothesis 1
p.adj <- c(p.adj,p.values[6]) # use unadjusted p-value for index

for (i in 1:length(outvars)){
  df <- data.frame(p=p.adj[i],estimate=estimates[i])
  plist3[[i]] <- plist3[[i]] +
    geom_text(data=df,aes(x=1,y=estimate,label=paste("p=",round(p,2))),col="black",size=3,nudge_y=0,nudge_x = 0.3) +
    geom_text(data=df,aes(x=1,y=estimate,label=paste("Diff=",round(estimate,2))),col="black",size=3,nudge_y=0.05,nudge_x = 0.3)
}

png(filename = "Figure_2.png")

gridExtra::grid.arrange(
  arrangeGrob(grobs = plist1,ncol=6,
              bottom=ggpubr::text_grob(sub("\n"," ",indvars.lab[1]),face="bold",vjust = -0.8, size=10)),
  arrangeGrob(grobs = plist2,ncol=6,
              bottom=ggpubr::text_grob(sub("\n"," ",indvars.lab[2]),face="bold",vjust = -0.8, size=10)),
  arrangeGrob(grobs = plist3,ncol=6,
              bottom=ggpubr::text_grob(sub("\n"," ",indvars.lab[3]),face="bold",vjust = -0.8, size=10)),
  left=textGrob("Difference in mean level of support",rot = 90, vjust = 1,gp = gpar(cex = 0.8))
)
dev.off()

# Figure 3: Factual manipulation checks --------

treat.rec.lab <- c("0 = otherwise, 1= Swedish citizens","0 = otherwise, 1= Women and children","0 = otherwise, 1= IS-affiliation")

# code adapted from Coppock 2020 replication kit
FMCplist1 <- list()
FMCplist2 <- list()

for (i in 1:length(FMCvars)){
  df <-
    PP[,c(indvars[i],FMCvars[i])] %>%
    gather(dv, value, FMCvars[i]) %>% mutate(IDV = factor(get(indvars[i]),
                                                          levels = c(0, 1),
                                                          labels = c("Otherwise", unlist(strsplit(indvars.lab[i],"\n"))[1])))
summary_df <-
    df %>%
    group_by(IDV, dv) %>%
    do(tidy(estimatr::lm_robust(value ~ 1, data = .))) %>%
    mutate(value = estimate, DV=FMCvars.lab[i])
  
FMCplist1[[i]] <- ggplot(summary_df, aes(IDV, value)) +
    geom_point(
      data = df,
      position = position_jitter(width = 0.2, height = 0.1),
      alpha = 0.03
    ) +
    geom_point(size = 1.5) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
    scale_y_continuous(breaks = seq(0, 1, length.out = 5)) +
    facet_wrap( ~ DV) +
    ylab(" ") +
    xlab(" ") 
  
diff_df <- broom::tidy(estimatr::lm_robust(get(FMCvars[i]) ~ get(indvars[i]), data = PP)) %>% mutate(IDV = indvars.lab[i], DV = FMCvars.lab[i]) %>% filter(term!="(Intercept)")
  
FMCplist2[[i]] <- ggplot(diff_df, aes(x = IDV, y = estimate)) + facet_wrap( ~ DV, scales = "fixed", shrink = TRUE, ncol = 3) +
    scale_y_continuous(limits=c(-0.0001,0.21)) +
    geom_point() +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
    geom_hline(yintercept = 0, colour = gray(0.5), lty = 2) + ylab("") + xlab("") 
  
}

png(filename = "Figure_3.png")

gridExtra::grid.arrange(
  arrangeGrob(grobs = FMCplist1,ncol=3,
              left=textGrob("1= received treatment\n0= otherwise",rot = 90, vjust = 1,gp = gpar(cex = 0.8)),
              bottom=ggpubr::text_grob("Randomly assigned identity factor",face="bold",vjust = -0.8, size=10)),
  arrangeGrob(grobs = FMCplist2,ncol=3,
              left=textGrob("Difference in probability",rot = 90, vjust = 1,gp = gpar(cex = 0.8)),
              bottom=ggpubr::text_grob("Comparison of randomly assigned identity factors",face="bold", size=10, vjust=-0.8)))

dev.off()




# Figure 4: Subjective manipulation checks --------------------------------

dflist <- list()

for (i in 1:length(SMCvars)){
  dflist[[i]] <- 
    rbind(broom::tidy(estimatr::lm_robust(get(SMCvars[i]) ~ SwedCitizens, data = PP)),
          broom::tidy(estimatr::lm_robust(get(SMCvars[i]) ~ WomenChildren, data = PP)), 
          broom::tidy(estimatr::lm_robust(get(SMCvars[i]) ~ ISaffiliation, data = PP))) %>% mutate(DV = SMCvars.lab[i])
}

df <- do.call(rbind,dflist) %>% mutate(DV = factor(DV, levels = SMCvars.lab)) %>% filter(term!="(Intercept)") %>% mutate(IDV = factor(term, levels = indvars,labels=indvars.lab))

png(filename = "Figure_4.png")

ggplot(df, aes(x = IDV, y = estimate)) + facet_wrap( ~ DV, scales = "fixed", shrink = TRUE, ncol = 2) + geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
  geom_hline(yintercept = 0, colour = gray(0.5), lty = 2) + ylab("Difference in probability") + xlab("Comparison of randomly assigned identity factors") 

dev.off()

# Figure 5: Placebo tests -------------------------------------------------

dflist <- list()

for (i in 1:length(PTvars)){
  dflist[[i]] <- rbind(broom::tidy(estimatr::lm_robust(get(PTvars[i]) ~ SwedCitizens, data = PP)),
                       broom::tidy(estimatr::lm_robust(get(PTvars[i]) ~ WomenChildren, data = PP)), 
                       broom::tidy(estimatr::lm_robust(get(PTvars[i]) ~ ISaffiliation, data = PP))) %>% mutate(DV = PTvars.lab[i])
}

df <- do.call(rbind,dflist) %>% mutate(DV = factor(DV, levels = PTvars.lab)) %>% filter(term!="(Intercept)") %>% mutate(IDV = factor(term, levels = indvars,labels=indvars.lab))

png(filename = "Figure_5.png")

ggplot(df, aes(x = IDV, y = estimate)) + facet_wrap( ~ DV, scales = "fixed", shrink = TRUE, ncol = 3) + geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
  geom_hline(yintercept = 0, colour = gray(0.5), lty = 2) + ylab("Difference in probability") + xlab("Comparison of randomly assigned identity factors") 

dev.off()

# Figure 6: Mediation analysis: moral obligation and threat perceptions --------

## mediator models
rm(df)

med.models <- list(update(covadj_model, Moral ~ + .),update(covadj_model, Threat ~ + .),update(covadj_model, Exp_Suc ~ + .),update(covadj_model, Exp_Cost ~ + .),update(covadj_model, Exp_Risk ~ + .),update(covadj_model, Exp_Ben ~ + .))

# omit missing for calculate of ACMEs
PP.nona <- na.omit(PP[,c(indvars,medvars,outvars,adjcovars)])


medlist1 <- list() # store lm mediator models for mediate (using na omit dataset)
medlist2 <- list() # one list for each IDV
medlist3 <- list() # each i is one mediating variable
meddf1 <- list() # dfs of robust lm mediator models for plotting
meddf2 <- list()  
meddf3 <- list() 

for (i in 1:length(medvars)){
  med.model <- med.models[[i]] # for some reason needs this
  medlist1[[i]] <- lm(update(med.model, . ~ SwedCitizens + .), data = PP.nona) 
  meddf1[[i]] <- broom::tidy(estimatr::lm_robust(update(med.model, . ~ SwedCitizens + .), data = PP)) %>% mutate(MV = medvars.lab[i])
  
  med.model <- med.models[[i]]
  medlist2[[i]] <- lm(update(med.model, . ~ WomenChildren + .), data = PP.nona)
  meddf2[[i]] <- broom::tidy(estimatr::lm_robust(update(med.model, . ~ WomenChildren + .), data = PP)) %>% mutate(MV = medvars.lab[i])
  
  med.model <- med.models[[i]]
  medlist3[[i]] <- lm(update(med.model, . ~ ISaffiliation + .), data = PP.nona)
  meddf3[[i]] <- broom::tidy(estimatr::lm_robust(update(med.model, . ~ ISaffiliation + .), data = PP)) %>% mutate(MV = medvars.lab[i])
}

# use this for summary plot
med_df <- rbind(do.call(rbind,meddf1),
                do.call(rbind,meddf2),
                do.call(rbind,meddf3)) %>% mutate(MV = factor(MV, levels = medvars.lab)) %>% filter(term %in% indvars) %>% mutate(IDV = factor(term, levels = indvars, labels=indvars.lab))


## outcome models
out.model <- update(covadj_model, Supp_Mil ~ . + Moral + Threat + Exp_Suc + Exp_Cost + Exp_Risk + Exp_Ben)

outlist1 <- list() # lm outcome models for mediate
outlist2 <- list() # one list for each IV
outlist3 <- list() # each i is one DV
outdf1 <- list()  
outdf2 <- list() 
outdf3 <- list()

for (i in 1:length(outvars)){
  outlist1[[i]] <- lm(update(out.model, get(outvars[i]) ~ SwedCitizens + . ), data = PP.nona)
  outdf1[[i]] <- broom::tidy(estimatr::lm_robust(update(out.model, get(outvars[i]) ~ SwedCitizens + . ), data = PP)) %>% mutate(IDV=indvars.lab[1]) %>% mutate(DV=outvars.lab[i])
  outlist2[[i]] <- lm(update(out.model, get(outvars[i]) ~ WomenChildren + . ), data = PP.nona)
  outdf2[[i]] <- broom::tidy(estimatr::lm_robust(update(out.model, get(outvars[i]) ~ WomenChildren + . ), data = PP)) %>% mutate(IDV=indvars.lab[2]) %>% mutate(DV=outvars.lab[i])
  outlist3[[i]] <- lm(update(out.model, get(outvars[i]) ~ ISaffiliation + . ), data = PP.nona)
  outdf3[[i]] <- broom::tidy(estimatr::lm_robust(update(out.model, get(outvars[i]) ~ ISaffiliation + . ), data = PP)) %>% mutate(IDV=indvars.lab[3]) %>% mutate(DV=outvars.lab[i])
}


# use this for summary plot
out_df1 <- rbind(do.call(rbind,outdf1),
                 do.call(rbind,outdf2),
                 do.call(rbind,outdf3)) %>% filter(term %in% medvars[1:2]) %>% mutate(MV = factor(term, levels = medvars[1:2],labels=medvars.lab[1:2])) %>% mutate(IDV = factor(IDV, levels = indvars.lab)) %>% mutate(DV = factor(DV, levels = outvars.lab))

## average causal mediation effects

med.model <- med.models[[1]] # for some reason totally necessary
med <- lm(update(med.model, . ~ SwedCitizens + .), data = PP.nona) 
out <- lm(update(out.model, . ~ SwedCitizens + .), data = PP.nona)
ace <- mediation::mediate(med,out,treat="SwedCitizens", mediator=medvars[1],sims=1000,robustSE = TRUE)
#summary(ace)

# prepare nested lists to store mediation models
ilist <- list("model1","model2")
jlist <- list("model3","model4","model5")
acelist1 <- list(ilist,jlist)
acelist2 <- list(ilist,jlist)
acelist3 <- list(ilist,jlist)

acedf1 <- list(ilist,jlist) # dfs with avg acmes and prop mediated added
acedf2 <- list(ilist,jlist)
acedf3 <- list(ilist,jlist)

# use only first two mediators (i)
# use all outcome variables (j)
# Use 1000 sims

for (i in 1:2){
  for (j in 1:length(outvars)){
    acelist1[[i]][[j]] <- mediation::mediate(medlist1[[i]],outlist1[[j]],treat="SwedCitizens", mediator=medvars[i],sims=100,robustSE = TRUE) 
    acedf1[[i]][[j]] <- broom::tidy(acelist1[[i]][[j]],conf.int = TRUE) %>% 
      add_row(term="acme",estimate=acelist1[[i]][[j]]$d.avg,std.error=0,p.value=acelist1[[i]][[j]]$d.avg.p,conf.low=acelist1[[i]][[j]]$d.avg.ci[1],conf.high=acelist1[[i]][[j]]$d.avg.ci[2]) %>% mutate(MV = medvars.lab[i]) %>% mutate(IDV = indvars.lab[1]) %>% mutate(DV = outvars.lab[j])
    
    acelist2[[i]][[j]] <- mediation::mediate(medlist2[[i]],outlist2[[j]],treat="WomenChildren", mediator=medvars[i],sims=100,robustSE = TRUE) 
    acedf2[[i]][[j]] <- broom::tidy(acelist2[[i]][[j]],conf.int = TRUE) %>% 
      add_row(term="acme",estimate=acelist2[[i]][[j]]$d.avg,std.error=0,p.value=acelist2[[i]][[j]]$d.avg.p,conf.low=acelist2[[i]][[j]]$d.avg.ci[1],conf.high=acelist2[[i]][[j]]$d.avg.ci[2]) %>% mutate(MV = medvars.lab[i]) %>% mutate(IDV = indvars.lab[2]) %>% mutate(DV = outvars.lab[j])
    
    acelist3[[i]][[j]] <- mediation::mediate(medlist3[[i]],outlist3[[j]],treat="ISaffiliation", mediator=medvars[i],sims=100,robustSE = TRUE) 
    acedf3[[i]][[j]] <- broom::tidy(acelist3[[i]][[j]],conf.int = TRUE) %>% 
      add_row(term="acme",estimate=acelist3[[i]][[j]]$d.avg,std.error=0,p.value=acelist3[[i]][[j]]$d.avg.p,conf.low=acelist3[[i]][[j]]$d.avg.ci[1],conf.high=acelist3[[i]][[j]]$d.avg.ci[2]) %>% mutate(MV = medvars.lab[i]) %>% mutate(IDV = indvars.lab[3]) %>% mutate(DV = outvars.lab[j])
  }}

acme_df1 <- rbind(do.call(rbind,acedf1[[1]]),
                  do.call(rbind,acedf1[[2]]),
                  do.call(rbind,acedf2[[1]]),
                  do.call(rbind,acedf2[[2]]),
                  do.call(rbind,acedf3[[1]]),
                  do.call(rbind,acedf3[[2]])) %>% mutate(MV = factor(MV, levels = medvars.lab[1:2])) %>% mutate(IDV = factor(IDV,levels=indvars.lab)) %>% filter(term=="acme") %>% mutate(DV = factor(DV,levels=outvars.lab))


# adjust p-values
acme_df1 <- acme_df1 %>% arrange(MV,DV) %>% mutate(p.adj = c(p.adjust(p.value[1:15],method="holm"),p.value[16:18],
                                                             p.adjust(p.value[19:33],method="holm"),p.value[34:36]))

# rename so that can be combined with other dfs for summary plot
acme_df1 <- acme_df1 %>% mutate(p.value= p.adj)


med_df <- rbind(med_df %>% mutate(DV=outvars.lab[1]),
                med_df %>% mutate(DV=outvars.lab[2]),
                med_df %>% mutate(DV=outvars.lab[3]),
                med_df %>% mutate(DV=outvars.lab[4]),
                med_df %>% mutate(DV=outvars.lab[5]),
                med_df %>% mutate(DV=outvars.lab[6])) %>% mutate(DV=factor(DV,levels=outvars.lab))

mediation1 <- rbind(med_df %>% dplyr::select(c("term","estimate","p.value","conf.low","conf.high","MV","IDV","DV")) %>% mutate(Model="Effect of identity factor\non mediator"),
                    out_df1 %>% dplyr::select(c("term","estimate","p.value","conf.low","conf.high","MV","IDV","DV")) %>% mutate(Model="Effect of mediator\non outcome"),
                    acme_df1 %>% dplyr::select(c("term","estimate","p.value","conf.low","conf.high","MV","IDV","DV"))  %>% mutate(Model="Average causal mediation effect")) %>% mutate(Model= factor(Model, levels= c("Effect of identity factor\non mediator","Effect of mediator\non outcome","Average causal mediation effect")))  %>% filter(MV %in% medvars.lab[1:2]) %>% mutate(MV = factor(MV,levels=medvars.lab[1:2],labels=paste("Mediator:",medvars.lab[1:2]))) %>% mutate(DV = factor(DV,levels=outvars.lab,labels=paste("Outcome: ",outvars.lab))) %>% mutate(sig=factor(ifelse(p.value<=0.05,1,0),levels=c(0,1)))

png(filename = "Figure_6.png")

ggplot(mediation1, aes(x = IDV, y = estimate,color=sig)) + facet_wrap(DV ~ MV, scales = "fixed", shrink = TRUE, ncol = 2) + geom_point(aes(shape=Model),position = position_dodge(width=0.5)) +
  scale_shape_manual(values=c(15:18)) +
  scale_color_manual(values = c("black", "red")) +
  geom_errorbar(aes(x=IDV,group=Model, ymin = conf.low, ymax = conf.high,color=sig), width = 0, position = position_dodge(width=0.5)) +
  geom_hline(yintercept = 0, colour = gray(0.5), lty = 2) + ylab("Estimate") + xlab("Comparison of randomly assigned identity factors") + guides(color="none")

dev.off()








